#!/usr/bin/perl -w
use strict;
use lib ("/net/cpp-group/Leo/bin");
use leo_GOcw qw{ mk_GO_annots_file_for_genes get_GO_slims read_GO_data analyse_GO_freqs};
use lib "/net/cpp-group/caleb/new_node_stuff/my_perl_stuff";
use cwstats_fgu202 qw{ do_hyperGlogcalc_this_or_greater do_hyperGlogcalc_this_or_fewer hyperGlogcalc};
use DBI;
use Getopt::Long;

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

my $usage = <<'USAGE';

USAGE:

    analyse_go_lite.pl
				--background FILE
				--sample FILE 
				--significance CUTOFF
				--temp_dir DIRECTORY
	analyses go terms for gene in --sample agains the --background
	prints to STDOUT
USAGE


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory

# optional parameters
my $background;
my $sample;
my $significance;
my $temp_dir;
my $help;
GetOptions(
			'background=s'			=> \$background,   
			'sample=s'				=> \$sample,       
			'significance=s'		=> \$significance, 
			'temp_dir=s'			=> \$temp_dir,     
			'help'					=> \$help);


die $usage if ($help);
die "Error:\nMissing parameter [$background].\n".$usage unless ($background);
die "Error:\nMissing parameter [$sample].\n".$usage unless ($sample);
die "Error:\nMissing parameter [$significance].\n".$usage unless ($significance);
die "Error:\nMissing parameter [$temp_dir].\n".$usage unless ($temp_dir);













my $host='db';
my $user='caleb';
my $database= 'ensembl_mart_26_1';
my $password='';
my $dbh_ens = DBI->connect("DBI:mysql:database=$database;host=$host",$user, $password,{'RaiseError' => 1});

#get background

my @bggenes;
open(BGSET,$background) || die "Error:\n\tCould not open [$background]\n$!";
while (<BGSET>)
{
	chomp;
	push @bggenes,$_;
}

my $numgenes = @bggenes;
print "numgenes $numgenes\n";

mk_GO_annots_file_for_genes_ensmart27_largegenesets(\@bggenes,$dbh_ens,"$temp_dir/temp_GO_lm_bg","$temp_dir/temp_GO_bp_bg","$temp_dir/temp_GO_fm_bg","hsapiens");
print "made files\n";

my (%goslim,%go_slim_all_ids);
get_GO_slims(\%goslim,\%go_slim_all_ids);

my (%gene2GO,%loc_freqs,%bioproc_freqs,%funcmap_freqs,%goid2annot);
read_GO_data(\%gene2GO,\%go_slim_all_ids,\%loc_freqs,\%bioproc_freqs,\%funcmap_freqs,\%goid2annot,"temp_GO_lm_bg","temp_GO_bp_bg","temp_GO_fm_bg");

open(SAMSET,$sample) || die "Error:\n\tCould not open [$sample]\n$!";
my @samgenes = <SAMSET>;
chomp @samgenes;







#_________________________________________________________________________________________
#
#
#	retrieve_GO_ids_per_category
#
#_________________________________________________________________________________________
sub retrieve_GO_ids($$$$\%\%)
{
	my ($dbh, $ens_species, $gene_id_to_GO_ids, $GO_id_to_description) = @_;
	
	my @category_data = 
	(
		["L" , "glook_biol_process_id", "$ens_species\_gene_ensembl__go_biol_process__look"]
		["BP", "glook_cell_location_id", "$ens_org\_gene_ensembl__go_cell_location__look"  ]
		["FM", "glook_mol_function_id", "$ens_org\_gene_ensembl__go_mol_function__look"    ]
	);
	for (@category_data)
	{
		my ($GO_cat, $field,  $table) = $_;
		retrieve_GO_ids_per_category($dbh, $table, $field, 
									 %{$gene_id_to_GO_ids->{$GO_cat}},
									 %$GO_id_to_description);
	}
	
}
#_________________________________________________________________________________________
#
#
#	retrieve_GO_ids_per_category
#
#_________________________________________________________________________________________
sub retrieve_GO_ids_per_category($$$\%\%)
{
	my ($dbh, $table_name, $field, $gene_id_to_GO_ids, $GO_id_to_description) = @_;
	my $sql_data = $dbh->selectall_arrayref(<<'PL/SQLCMD');
        SELECT 
            gene_stable_id AS gene_id, 
            $field AS GO_id, 
            description
        FROM 
            $table_name
        WHERE
            NOT $field is null;
PL/SQLCMD
	for my $line(@$sql_data)
	{
		push(@{$gene_id_to_GO_ids->{$line->[0]}}, $line->[1]);
		$GO_id_to_description->{$line->[1]} = $line->[2];
	}

}






#
#	print GO terms which are under or over represented
#
sub print_GO_over_under($$\@\@\@)
{
	my ($dbh, $ens_species, $gene_ids, $gene_ids_background, $significant_results_all) = @_;

	# get GO identifiers and descriptions
	my %gene_id_to_GO_ids_all;
	my %GO_id_to_description;
	retrieve_GO_ids($dbh, $ens_species, %gene_id_to_GO_ids_all, %GO_id_to_description);
	
	# get GO slim identifiers
	my %GO_slim_ids;
	get_GO_slims(%GO_slim_ids, "/net/cpp-group/caleb/new_node_stuff/GO/goslim_generic.obo");

	# get significant results
	my %significant_results;
	get_under_overrepresented_GO_terms(	%gene_id_to_GO_ids_all, %GO_slim_ids, 
										 @$gene_ids_background, @$gene_ids, 
										 0.01, %significant_results);


	# print out for each GO category
	for my $GO_cat (keys %significant_results)
	{
		# for each set of under and then over represented
		for my $over_under (keys %{$significant_results->{$GO_cat}})
		{
			# sort GO identifiers by their descriptions
			my @GO_ids = sort {	$GO_id_to_description{$a} cmp
								$GO_id_to_description{$b}}
								keys %{$significant_results->{$GO_cat}{$over_under}};

			# for each GO identifier
			for my $GO_id(@GO_ids)
			{
				printf	"$GO_cat\t$over_under\t$GO_id\t%d\t%d\t%d\t%d\t%1.2e\t%1.2e\n",
							@{$significant_results->{$GO_cat}{$over_under}{$GO_id}};
			}
		}
		
	}
}

